% By Martin Andreasen, 2009
% This function computes bond prices up to a given maturity, provided the
% first bond prices is reported in the function g(x,sig) and stored at the
% position ny_1bond. This implementation allows for no-zero third moment of
% the shock distribution.
% The integer "R_function" has the following properties:
% 1  - if log-transformation of bond prices is used and 
% 0  - if no log-transformation of bond prices is used

% In this version, symmetries in Pxx and Pxxx are used

function [P,Px,Pxx,Pxxx,Pss,Pssx,Psss] = Get_Bond_Prices_Mom3(M,Mxp,Myp,...
    Mxpx,Mxpxp,Mxpy,Mxpyp,Mypx,Mypxp,Mypy,Mypyp,...
    gx,gxx,gxxx,gss,gssx,gsss,hx,hxx,hxxx,hss,hssx,hsss,eta,vectorMom3,Maturity,ny_1bond,R_function,order_app);


% The arrays
[ny,nx] = size(gx);
ne      = size(eta,2);
P       = zeros(Maturity,1);
Px      = zeros(Maturity,nx);
Pxx     = zeros(Maturity,nx,nx);
Pxxx    = zeros(Maturity,nx,nx,nx);
Pss     = zeros(Maturity,1);
Pssx    = zeros(Maturity,nx);
Psss    = zeros(Maturity,1);

% The all the derivatives for the first bond price
P(1,1)       = M;
Px(1,:)      = gx(ny_1bond,:);
Pxx(1,:,:)   = squeeze(gxx(ny_1bond,:,:));
Pxxx(1,:,:,:)= squeeze(gxxx(ny_1bond,:,:,:));
Pss(1,1)     = gss(ny_1bond,1);
Pssx(1,:)    = gssx(ny_1bond,:);
Psss(1,1)    = gsss(ny_1bond,1);

% The value of bond prices in deterministic steady state steady state
k = 2:Maturity;
P(k,1) = M.^k;
%for k=2:Maturity;
%    P(k,1) = M^k;
%end

if R_function == 2
    % The logit transformation is used for bonds: p to 1/(exp(phat)+1)
    PhatSS    = log(1./P(:,1)-1);
    R(:,1)    = (exp(PhatSS(:,1))+1).^(-1);
    Rp(:,1)   = -(exp(PhatSS(:,1))+1).^(-2).*exp(PhatSS(:,1));
    Rpp(:,1)  = 2*(exp(PhatSS(:,1))+1).^(-3).*exp(2*PhatSS(:,1))+Rp(:,1);
    % ERROR here. Check if formula for Pxxx is correct!!!!
    Rppp(:,1) = -6*(exp(PhatSS(:,1))+1).^(-4).*exp(3*PhatSS(:,1))+4*(exp(PhatSS(:,1))+1).^(-3).*exp(2*PhatSS(:,1))+Rpp(:,1);
elseif R_function == 1
    % Log-transformation is used of bond prices: p to exp(phat)
    R(:,1)    = P(:,1);
    Rp(:,1)   = R(:,1);
    Rpp(:,1)  = R(:,1);    
    Rppp(:,1) = R(:,1);    
elseif R_function == 0
    % No transformation is used of bond prices: p to phat    
    R(:,1)   = P(:,1);
    Rp(:,1)  = ones(Maturity,1);
    Rpp(:,1) = zeros(Maturity,1);    
    Rppp(:,1)= zeros(Maturity,1);     
end

% The first order terms
% Recall that R(P(0,1)) = 1 
for k=2:Maturity
   % With matrix notation
   Px(k,:) = (Px(1,:)*Rp(1,1))*R(k-1,1)+R(1,1)*Rp(k-1,1)*Px(k-1,:)*hx;
   Px(k,:) = Px(k,:)/Rp(k,1);
end
% Stop if order_app == 1
if order_app == 1
    return;
end

% The second order terms, Pxx
for k=2:Maturity
    for alfa1=1:nx
        for alfa2=alfa1:nx
            Pxx(k,alfa1,alfa2) = -Rpp(k,1)*Px(k,alfa2)*Px(k,alfa1) ... 
                +(Rp(1,1)*Pxx(1,alfa1,alfa2)+Rpp(1,1)*Px(1,alfa2)*Px(1,alfa1))*R(k-1,1)...
                +Px(1,alfa1)*Rp(1,1)*Rp(k-1,1)*Px(k-1,:)*hx(:,alfa2)...
                +Px(1,alfa2)*Rp(1,1)*Rp(k-1,1)*Px(k-1,:)*hx(:,alfa1)...
                +R(1,1)*Rpp(k-1,1)*Px(k-1,:)*hx(:,alfa2)*Px(k-1,:)*hx(:,alfa1)...
                +R(1,1)*Rp(k-1,1)*hx(:,alfa1)'*squeeze(Pxx(k-1,:,:))*hx(:,alfa2)...
                +R(1,1)*Rp(k-1,1)*Px(k-1,:)*squeeze(hxx(:,alfa1,alfa2));
        end
        if alfa1 > 1
            Pxx(k,alfa1,1:alfa1-1) = squeeze(Pxx(k,1:alfa1-1,alfa1))';
        end
    end
    Pxx(k,:,:) = Pxx(k,:,:)/Rp(k,1);
end
  
% The constant correction, Pss
I_ne     = eye(ne);
for k=2:Maturity
   tmp = 0;
   for phi1=1:ne
      tmp = tmp + 2*Myp(1,:)*gx(:,:)*eta(:,phi1)*Rp(k-1,1)*Px(k-1,:)*eta(:,:)*I_ne(:,phi1)...
                + 2*Mxp(1,:)*eta(:,phi1)*Rp(k-1,1)*Px(k-1,:)*eta(:,:)*I_ne(:,phi1)...
                + R(1,1)*Rpp(k-1,1)*Px(k-1,:)*eta(:,phi1)*Px(k-1,:)*eta(:,:)*I_ne(:,phi1)...
                + R(1,1)*Rp(k-1,1)*eta(:,phi1)'*squeeze(Pxx(k-1,:,:))*eta(:,:)*I_ne(:,phi1);
   end
   Pss(k,1) = Pss(1,1)*Rp(1,1)*R(k-1,1) + tmp...
        + R(1,1)*Rp(k-1,1)*Px(k-1,:)*hss(:,1)...
        + R(1,1)*Rp(k-1,1)*Pss(k-1,1);

    Pss(k,1) = Pss(k,1)/Rp(k,1);
end
% Stop if order_app == 2
if order_app == 2
    return;
end


% The third order order terms, Pxxx
for k=2:Maturity
    for alfa1=1:nx
        for alfa2=alfa1:nx
            for alfa3=alfa2:nx
                tmp = 0;
                for gama1=1:nx
                    tmp = tmp + hx(gama1,alfa1)*hx(:,alfa2)'*(squeeze(Pxxx(k-1,gama1,:,:)))*hx(:,alfa3);
                end
                Pxxx(k,alfa1,alfa2,alfa3) = -Rppp(k,1)*Px(k,alfa3)*Px(k,alfa2)*Px(k,alfa1)...                           %1)
                    -Rpp(k,1)*Pxx(k,alfa2,alfa3)*Px(k,alfa1)...                                                         %2)
                    -Rpp(k,1)*Px(k,alfa2)*Pxx(k,alfa1,alfa3)...                                                         %3)
                    -Rpp(k,1)*Px(k,alfa3)*Pxx(k,alfa1,alfa2)...                                                         %4)
                    +(Rp(1,1)*Pxxx(1,alfa1,alfa2,alfa3)+Rppp(1,1)*Px(1,alfa3)*Px(1,alfa2)*Px(1,alfa1)+Rpp(1,1)*Pxx(1,alfa2,alfa3)*Px(1,alfa1)... %5)
                      +Rpp(1,1)*Px(1,alfa2)*Pxx(1,alfa1,alfa3)+Rpp(1,1)*Px(1,alfa3)*Pxx(1,alfa1,alfa2))*R(k-1,1)...     %6)
                    +(Rp(1,1)*Pxx(1,alfa1,alfa2)+Rpp(1,1)*Px(1,alfa2)*Px(1,alfa1))*Rp(k-1,1)*Px(k-1,:)*hx(:,alfa3)...   %7)
                    +(Rp(1,1)*Pxx(1,alfa1,alfa3)+Rpp(1,1)*Px(1,alfa3)*Px(1,alfa1))*Rp(k-1,1)*Px(k-1,:)*hx(:,alfa2)...   %8)
                    +Px(1,alfa1)*Rp(1,1)*Rpp(k-1,1)*Px(k-1,:)*hx(:,alfa3)*Px(k-1,:)*hx(:,alfa2)...                      %9)
                    +Px(1,alfa1)*Rp(1,1)*Rp(k-1,1)*hx(:,alfa2)'*squeeze(Pxx(k-1,:,:))*hx(:,alfa3)...                    %10)
                    +Px(1,alfa1)*Rp(1,1)*Rp(k-1,1)*Px(k-1,:)*squeeze(hxx(:,alfa2,alfa3))...                             %11)
                    +(Rp(1,1)*Pxx(1,alfa2,alfa3)+Rpp(1,1)*Px(1,alfa3)*Px(1,alfa2))*Rp(k-1,1)*Px(k-1,:)*hx(:,alfa1)...   %12)
                    +Px(1,alfa2)*Rp(1,1)*Rpp(k-1,1)*Px(k-1,:)*hx(:,alfa3)*Px(k-1,:)*hx(:,alfa1)...                      %13)
                    +Px(1,alfa2)*Rp(1,1)*Rp(k-1,1)*hx(:,alfa1)'*squeeze(Pxx(k-1,:,:))*hx(:,alfa3)...                    %14)
                    +Px(1,alfa2)*Rp(1,1)*Rp(k-1,1)*Px(k-1,:)*squeeze(hxx(:,alfa1,alfa3))...                             %15)
                    +Px(1,alfa3)*Rp(1,1)*Rpp(k-1,1)*Px(k-1,:)*hx(:,alfa2)*Px(k-1,:)*hx(:,alfa1)...                      %16)
                    +R(1,1)*Rppp(k-1,1)*Px(k-1,:)*hx(:,alfa3)*Px(k-1,:)*hx(:,alfa2)*Px(k-1,:)*hx(:,alfa1)...            %17)
                    +R(1,1)*Rpp(k-1,1)*hx(:,alfa2)'*squeeze(Pxx(k-1,:,:))*hx(:,alfa3)*Px(k-1,:)*hx(:,alfa1)...          %18)
                    +R(1,1)*Rpp(k-1,1)*Px(k-1,:)*squeeze(hxx(:,alfa2,alfa3))*Px(k-1,:)*hx(:,alfa1)...                   %19)
                    +R(1,1)*Rpp(k-1,1)*Px(k-1,:)*hx(:,alfa2)*hx(:,alfa1)'*squeeze(Pxx(k-1,:,:))*hx(:,alfa3)....         %20)
                    +R(1,1)*Rpp(k-1,1)*Px(k-1,:)*hx(:,alfa2)*Px(k-1,:)*hxx(:,alfa1,alfa3)...                            %21)
                    +Px(1,alfa3)*Rp(1,1)*Rp(k-1,1)*hx(:,alfa1)'*squeeze(Pxx(k-1,:,:))*hx(:,alfa2)...                    %22)
                    +R(1,1)*Rpp(k-1,1)*Px(k-1,:)*hx(:,alfa3)*hx(:,alfa1)'*squeeze(Pxx(k-1,:,:))*hx(:,alfa2)...          %23)
                    +R(1,1)*Rp(k-1,1)*tmp...                                                                            %24)
                    +R(1,1)*Rp(k-1,1)*hx(:,alfa1)'*squeeze(Pxx(k-1,:,:))*squeeze(hxx(:,alfa2,alfa3))....                %25)
                    +R(1,1)*Rp(k-1,1)*reshape(hxx(:,alfa1,alfa3),1,nx)*squeeze(Pxx(k-1,:,:))*hx(:,alfa2)....            %26)
                    +Px(1,alfa3)*Rp(1,1)*Rp(k-1,1)*Px(k-1,:)*squeeze(hxx(:,alfa1,alfa2))...                             %27)
                    +R(1,1)*Rpp(k-1,1)*Px(k-1,:)*hx(:,alfa3)*Px(k-1,:)*squeeze(hxx(:,alfa1,alfa2))...                   %28)
                    +R(1,1)*Rp(k-1,1)*reshape(hxx(:,alfa1,alfa2),1,nx)*squeeze(Pxx(k-1,:,:))*hx(:,alfa3)...             %29)
                    +R(1,1)*Rp(k-1,1)*Px(k-1,:)*squeeze(hxxx(:,alfa1,alfa2,alfa3));                                     %30)
                
                % Using symmetry for alfa1 and alfa2
                if alfa1 == alfa2 && alfa2 ~= alfa3 %alfa1==alfa2~=alfa3
                    %Pxxx(k,alfa1,alfa1,alfa3)= Pxxx(k,alfa1,alfa2,alfa3);
                    Pxxx(k,alfa1,alfa3,alfa1) = Pxxx(k,alfa1,alfa2,alfa3);
                    Pxxx(k,alfa3,alfa1,alfa1) = Pxxx(k,alfa1,alfa2,alfa3);                
                end
                % Using symmetry for alfa2 and alfa3            
                if alfa1 ~= alfa2 && alfa2 == alfa3  %alfa1~=alfa2==alfa3 
                    %Pxxx(k,alfa1,alfa2,alfa2)= Pxxx(k,alfa1,alfa2,alfa3);                
                    Pxxx(k,alfa2,alfa1,alfa2) = Pxxx(k,alfa1,alfa2,alfa3);                
                    Pxxx(k,alfa2,alfa2,alfa1) = Pxxx(k,alfa1,alfa2,alfa3);                                
                end            
                % Using symmetry for alfa1,alfa2, and alfa3            
                if alfa1 ~= alfa2 && alfa1 ~= alfa3 &&  alfa2 ~= alfa3 %alfa1~=alfa2~=alfa3
                    %Pxxx(k,alfa1,alfa2,alfa3) = Pxxx(k,alfa1,alfa2,alfa3);
                    Pxxx(k,alfa1,alfa3,alfa2) = Pxxx(k,alfa1,alfa2,alfa3);
                    Pxxx(k,alfa3,alfa1,alfa2) = Pxxx(k,alfa1,alfa2,alfa3);
                    Pxxx(k,alfa3,alfa2,alfa1) = Pxxx(k,alfa1,alfa2,alfa3);
                    Pxxx(k,alfa2,alfa3,alfa1) = Pxxx(k,alfa1,alfa2,alfa3); 
                    Pxxx(k,alfa2,alfa1,alfa3) = Pxxx(k,alfa1,alfa2,alfa3);                 
                end
            end
        end
    end
    Pxxx(k,:,:,:) = Pxxx(k,:,:,:)/Rp(k,1);
end

% The third order order terms, Pssx
% First some auxiliary variables
Myp_gxx = zeros(nx,nx);
for gama1=1:nx
    for gama2=1:nx
        Myp_gxx(gama1,gama2) = Myp(1,:)*squeeze(gxx(:,gama1,gama2));
    end
end

for k=2:Maturity
    for alfa3=1:nx
        Pssx(k,alfa3) = -Rpp(k,1)*Px(k,alfa3)*Pss(k,1)...                                                                 %1)
            +(Rp(1,1)*Pssx(1,alfa3)+Rpp(1,1)*Px(1,alfa3)*Pss(1,1))*R(k-1,1)...                                            %2)
            +Pss(1,1)*Rp(1,1)*Rp(k-1,1)*Px(k-1,:)*hx(:,alfa3);                                                            %3)
        
        % Terms 4) to 12)
        vec_dyp(1:ny,1) = (Mypyp(:,:)*gx*hx(:,alfa3)+Mypy(:,:)*gx(:,alfa3)+Mypxp(:,:)*hx(:,alfa3)+Mypx(:,alfa3));     %4')
        vec_dxp(1:nx,1) = (Mxpyp(:,:)*gx*hx(:,alfa3)+Mxpy(:,:)*gx(:,alfa3)+Mxpxp(:,:)*hx(:,alfa3)+Mxpx(:,alfa3));     %7')        
        tmp = 0;
        for phi1=1:ne
            tmp = tmp + 2*vec_dyp'*gx*eta(:,phi1)*Rp(k-1,1)*Px(k-1,:)*eta*I_ne(:,phi1)...                                 %4) and 5)
                      + 2*eta(:,phi1)'*Myp_gxx(:,:)*hx(:,alfa3)*Rp(k-1,1)*Px(k-1,:)*eta*I_ne(:,phi1)...                  %6)
                      + 2*vec_dxp'*eta(:,phi1)*Rp(k-1,1)*Px(k-1,:)*eta*I_ne(:,phi1)...                                    %7) and 8)
                      + 2*Myp(1,:)*gx*eta(:,phi1)*Rpp(k-1,1)*Px(k-1,:)*hx(:,alfa3)*Px(k-1,:)*eta*I_ne(:,phi1)...         %9)
                      + 2*Mxp(1,:)*eta(:,phi1)*Rpp(k-1,1)*Px(k-1,:)*hx(:,alfa3)*Px(k-1,:)*eta*I_ne(:,phi1)...            %10)
                      + 2*Myp(1,:)*gx*eta(:,phi1)*Rp(k-1,1)*(eta(:,:)*I_ne(:,phi1))'*squeeze(Pxx(k-1,:,:))*hx(:,alfa3)...%11)
                      + 2*Mxp(1,:)*eta(:,phi1)*Rp(k-1,1)*(eta(:,:)*I_ne(:,phi1))'*squeeze(Pxx(k-1,:,:))*hx(:,alfa3);     %12)
        end
        Pssx(k,alfa3) = Pssx(k,alfa3) + tmp;
        
        
        % Term 13) to 18)
        tmp = 0;
        for phi2=1:ne
           tmp = tmp + Px(1,alfa3)*Rp(1,1)*Rpp(k-1,1)*Px(k-1,:)*eta(:,phi2)*Px(k-1,:)*eta*I_ne(:,phi2)...                 %13)
                     + R(1,1)*Rppp(k-1,1)*Px(k-1,:)*hx(:,alfa3)*Px(k-1,:)*eta(:,phi2)*Px(k-1,:)*eta*I_ne(:,phi2)...       %14)
                     + R(1,1)*Rpp(k-1,1)*eta(:,phi2)'*squeeze(Pxx(k-1,:,:))*hx(:,alfa3)*Px(k-1,:)*eta*I_ne(:,phi2)...     %15)
                     + R(1,1)*Rpp(k-1,1)*Px(k-1,:)*eta(:,phi2)*(eta*I_ne(:,phi2))'*squeeze(Pxx(k-1,:,:))*hx(:,alfa3)...   %16)
                     + Px(1,alfa3)*Rp(1,1)*Rp(k-1,1)*(eta*I_ne(:,phi2))'*squeeze(Pxx(k-1,:,:))*eta(:,phi2)...             %17)
                     + R(1,1)*Rpp(k-1,1)*Px(k-1,:)*hx(:,alfa3)*(eta*I_ne(:,phi2))'*squeeze(Pxx(k-1,:,:))*eta(:,phi2);     %18)
        end
        Pssx(k,alfa3) = Pssx(k,alfa3) + tmp;
        
        % Term 19)
        tmp = 0;
        for gama1=1:nx
            for phi2=1:ne
                tmp = tmp + R(1,1)*Rp(k-1,1)*eta(:,phi2)'*squeeze(Pxxx(k-1,gama1,:,:))*hx(:,alfa3)*eta(gama1,:)*I_ne(:,phi2);
            end
        end
        Pssx(k,alfa3) = Pssx(k,alfa3) + tmp;                
        
        % Term 20) to 25)
        tmp = 0;
        tmp = tmp + Px(1,alfa3)*Rp(1,1)*Rp(k-1,1)*Px(k-1,:)*hss(:,1)...                                                     %20)
                  + R(1,1)*Rpp(k-1,1)*Px(k-1,:)*hx(:,alfa3)*Px(k-1,:)*hss(:,1)...                                           %21)
                  + R(1,1)*Rp(k-1,1)*hss(:,1)'*squeeze(Pxx(k-1,:,:))*hx(:,alfa3)...                                         %22)
                  + R(1,1)*Rp(k-1,1)*Px(k-1,:)*hssx(:,alfa3)...                                                             %23)
                  + Px(1,alfa3)*Rp(1,1)*Rp(k-1,1)*Pss(k-1,1)...                                                             %24)
                  + R(1,1)*Rpp(k-1,1)*Px(k-1,:)*hx(:,alfa3)*Pss(k-1,1)...                                                   %25)
                  + R(1,1)*Rp(k-1,1)*Pssx(k-1,:)*hx(:,alfa3);                                                               %26)
        Pssx(k,alfa3) = Pssx(k,alfa3) + tmp;                              
    end
    Pssx(k,:) = Pssx(k,:)/Rp(k,1); 
end

% The third order order terms, Psss
for k=2:Maturity
    % Term b1)
    Psss(k,1) = Rp(1,1)*Psss(1,1)*R(k-1,1);                                                                                 %b1)

    % The terms b2), b3), b5), b6), b7), b8) and b9)
    tmp = 0;
    for phi1=1:ne
        tmp = tmp + 3*(gx*eta(:,phi1)*vectorMom3(phi1))'*Mypyp*gx*eta(:,phi1)*Rp(k-1,1)*Px(k-1,:)*eta(:,phi1)...            %b2)
                  + 6*(gx*eta(:,phi1)*vectorMom3(phi1))'*Mypxp*eta(:,phi1)*Rp(k-1,1)*Px(k-1,:)*eta(:,phi1)...               %b3)
                  + 3*(eta(:,phi1)*vectorMom3(phi1))'*Mxpxp*eta(:,phi1)*Rp(k-1,1)*Px(k-1,:)*eta(:,phi1)...                  %b5)
                  + 3*(Myp(1,:)*gx*eta(:,phi1) + Mxp(1,:)*eta(:,phi1))*vectorMom3(phi1)*Rpp(k-1,1)*Px(k-1,:)*eta(:,phi1)... %b6)
                     *Px(k-1,:)*eta(:,phi1)...
                  + 3*(Myp(1,:)*gx*eta(:,phi1) + Mxp(1,:)*eta(:,phi1))*vectorMom3(phi1)*Rp(k-1,1)...                        %b7)
                     *eta(:,phi1)'*squeeze(Pxx(k-1,:,:))*eta(:,phi1)...
                  + R(1,1)*Rppp(k-1,1)*Px(k-1,:)*eta(:,phi1)*Px(k-1,:)*eta(:,phi1)...                                       %b8)
                     *Px(k-1,:)*eta(:,phi1)*vectorMom3(phi1)...
                  + 3*R(1,1)*Rpp(k-1,1)*Px(k-1,:)*eta(:,phi1)...                                                            %b9)
                     *(eta(:,phi1)*vectorMom3(phi1))'*squeeze(Pxx(k-1,:,:))*eta(:,phi1);                           
    end
    Psss(k,1) = Psss(k,1) + tmp;
    
    % Term b4)
    tmp = 0;
    for beta1=1:ny
        for phi1=1:ne
             tmp = tmp + 3*(eta(:,phi1)*vectorMom3(phi1))'*Myp(1,beta1)*squeeze(gxx(beta1,:,:))*eta(:,phi1)...
                             *Rp(k-1,1)*Px(k-1,:)*eta(:,phi1);
        end
    end
    Psss(k,1) = Psss(k,1) + tmp;
    
    % The term 10b)
    tmp = 0;
    for gama1=1:nx
        for phi1=1:ne
            tmp = tmp + R(1,1)*Rp(k-1,1)*eta(:,phi1)'*squeeze(Pxxx(k-1,gama1,:,:))*eta(:,phi1)*eta(gama1,phi1)*vectorMom3(phi1);
        end
    end
    Psss(k,1) = Psss(k,1) + tmp;
 
    % Terms b11 and b12
    tmp =  R(1,1)*Rp(k-1,1)*Px(k-1,:)*hsss(:,1)...                                                              %b11)
         + R(1,1)*Rp(k-1,1)*Psss(k-1,1);                                                                        %b12)
    Psss(k,1) = Psss(k,1) + tmp;

    % Scaling by Rp
    Psss(k,1) = Psss(k,1)/Rp(k,1);
end

